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This work is the result of three years' effort at 
University of Virginia to assist NASA in choosing the numer- 
ical integrators to be used in real time simulators for 
aircraft and spacecraft. The first author was introduced to 
this subject during a visit to the Institute for Computer 
Applications in Science and Engineering at the NASA Langley 
Research Center during Summer 1975. This was followed by a 
grant NSG-1335 from 1976 to 1979. This report described 
three interrelated software systems written under this 
grant. 

The principal investigator wishes to thank his co- 
authors/graduate students, also students Sandra Bollinger 
and Bob Athay who did some of the initial work. Dr. 

R. L. Bowles of NASA both got us started and provided tech- 
nical guidance along the way. 

Copies of the programs for CDC equipment and IBM 360 or 
370 are available from the first author. For other machines, 
the IBM version is almost ANSI standard. 
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ABSTRACT 


Three computer programs for the evaluation and testing 
of numerical integration formulas for use with fixed step- 
size programs to solve initial value systems of ordinary 
differential equations are described. SERIES, written in 
PASCAL, takes as input the differential equations and pro- 
duces FORTRAN subroutines for the derivatives of the system 
and for computing the actual solution through recursive 
power series techniques. Both of these are used by STAN, a 
FORTRAN program to interactively display a discrete analog 
of the Liapunov stability region of any two-dimensional 
subspace of the system. The derivatives may be used by 
CLMP, a FORTRAN program, to test the fixed stepsize formula 
against a good numerical result and interactively display 
the solutions. 
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1 . INTRODUCTION 


Consider the numerical solution of the problem 

y (t) = f(y,t) 

y(0) = y 0 (i) 

for f:R n+1 •* R n , where f is a vector valued function of 
independent variable t and the dependent variable y(t) in 
R n , y Q is the initial value from R n , and y(t) is often 
called the state vector. This form can be used to study a 
wide variety of initial value problems including the fol- 
lowing: the Method of Lines (MOL) appioach to solving ini- 
tial-boundary value problems in parabolic systems of partial 
differential equations [Ames]; higher order ordinary dif- 
ferential equations using a change of variable 

y i+ i = d 1 y/dt 1 

3 3 

so that y' 3 = d y/dt [Gear]; mixed differential and non- 
linear algebraic systems F(y,y',t) = 0, where F (y,y',t) is 
solved by an implicit numerical method using Newton's iter- 
ative method or similar [Gear]; and the equations of state, 
including conservation laws, of an engineering simulation. 
A further simplification could replace the independent vari- 
able t by a new element y n+ ^ such that y' n+1 = 1 , y n+1 (0) = 
0. However, this will not be used in the present analysis. 

If (1) were a linear homogeneous equation 



then the solution would be trivial, as the following anal- 
ysis shows. For a matrix A of full rank there exists an 
Hermitian transformation P such that: 

PAP* = D 

where D is a diagonal matrix X^] of the eigenvalues of 
A. Write (2) as 

Py» = PA p* Py (3) 

Set z ( t ) * Py(t) so that z’ = Dz is equivalent to (2) and 
remember that z is complex since the eigenvalues of A may be 
complex. Then the solution of (2) is obtained from (3) as 

z i (t) = exp (\ i t)z i (0), (3) 

y(t) = P* z(t). 

It has been customary to investigate the stability and 
accuracy of numerical solutions of (1) by describing the 
effect of solving (2) numerically, or more simply, by in- 
vestigating the complex equation 

y = Ay, 

y(0) = y 0 (4) 

since (2) can be reduced to (4) in each component. This 
linear stability analysis will be described here, and its 
shortcomings pointed out, as a prelude to the description in 
Chapter 2 of the mathematical foundations of software speci- 
fically designed to analyze and test numerical methods on 
small (< 20 dependent variables) subsystems of nonlinear 
differential equations. 
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Common numerical formulae for integrating the system 
(1) are listed in Appendix A; they are either 1-step methods 

y n a ♦ <y n -r f ' v h) 

or multistep methods 

0 = L(y n' y n-l'**" y n-k' f » V 
where y n = y(0+nh) + e n is the approximate solution after n 

equal steps of size h in the independent variable t from the 

initial point ( y q, 0). The global error e n> which depends 

on h, the function f(y,t), and the numerical method is 

described elsewhere [Gear, Stetter]. However, the stability 

analysis is reviewed here. 

One-step methods are typified by explicit Runge-Kutta 
formulas 

K 0 = hf < y n-l' W' 

K q = hf < y n-l + f Q b qj K j' < 5 > 

q = 1, . . .,s-l, 

6-1 

y = y n , + I g„ K . 
n n-1 q=o <3 <3 

which attempt to approximate the Taylor series of y(t) about 
t n-1 by using a linear combination of S stages. 

Typical examples include the two-stage formula 

K 0 = hf (y n-l' t n-l ) 

K 1 s hf < y n+l + * 5 * K 0' Vl + h/2) 

y n = ^n-l + K 1- 
3 



This and the typical four-stage Runge-Kutta formula are 
given in Appendix A. 

The result in the linear case is that 
y n = (1 ♦ IA + (hA> 2 /2 + — (hXjP/pMv^j + 0(h p+1 ) (6) 

and the formula is called of error order p<s. Even if b g j 
and g q are picked so that p is as large as possible, the 
coefficients are still not fully specified and an infinite 
family of coefficients, in one or more parameters, results. 
In the linear case, all choices of (5) result in (6), so the 
linear stability analysis can be simply stated: 
let e n =■ y n - z n where y Q f z Q are initial values used in 
the numerical solution, then the stability region S = (hA: 
lim (n+») e n is finite). Since e(t n ) = exp (nhAt) (y Q - 
z Q ), S will contain values of hA where h>0 and Re (A) < 0 in 
the exact case. For (6), S will contain hA such that 
|l + hA +— +(hA) p /p! | < 1. 

Figure 1.1 shows the stability region of Runge-Kutta for- 
mulas of various order, and [Jeltsch] has shown that even if 
the coefficients are not picked to maximize p but rather to 
maximize the radius of the largest circle in S tangent to 
the imaginary axis at 0., this region is still bounded for 
finite s. 

In the multistep case, predictor only and piidictor- 
corrector combinations are used. Predictor formulas have 
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Figure 1.1 Stability regions of Runge-Kutta methods of 
order p«l ,2,3,5, (increasing area and 
7 (dotted outline). 
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Figure 1.2. Stability regions for Adams-Bashforth 
methods. Method of order K is stable 
inside region to the left of origin. 
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the form 


*n s y n -j + “j f <y n -j' Vj>) 


(7) 


after applying an appropriate predictor, an implicit correc- 


tor formula may be solved to improve the solution y R : 

0 * ilo y "-i + “l Vj> 

This can be rewritten (assuming a^ - -1.) as 

y n = 8 + ^ f <y«' tn) (?) 


n' n' 


where 


( 8 ) 


s = 


j=l (a 3 y n-3 + f <y n-j' 


(9) 


doesn't change during the iteration. One can solve (9) 

iteratively by letting = yj!| and solving 

y^) = s + hb!I f(y^ i ’ 1 ^, t ) 

*n 0 vjr n ' n' 

for l = 1, 2, . . . ,Q where Q is either set to a constant 
(often 1), or chosen after some convergence criterion such 
as (y^*- y^ i-1 > | < some small value. Newton's method can 
be used on the function 

° = y n - • - hbj f(y n , t n ) 

by iterating 

- <I-l*Jj n r 1 <y< i ‘ 1, -*-hbJ V* (10) 

where J n is a clotc approximation to the Jacobian matrix 
3f/3yf 


y.t- . 
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The linear stability analysis often assumes that (9) is 
solved exactly, but in fact if Q, the number of corrector 
iterations, is finite, this will affect the stability. If 
Newton's method is used with the exact Jacobian X, then only 
one iteration will yield the exact solution of (9), so this 
is valid in the linear case. Note that (8) is the linear, 
homogeneous difference equation 

d- b o »>r B y : ( *j ♦ b j M )y n -j 

and, for each particular hX, the closed- form solution can be 
obtained [Gear] . This solution is of the form 




k ra * j ^ 

I ( 2 1 w. . n 3 " 1 ) x? 
i=l j=l 13 1 


( 11 ) 


where x^ is one of the S < k unique roots of the polynomial 
equation 

k * < 

0 = I (a- + hXb*) x 3 , 
j=0 3 3 

ra^ is the multiplicity of the i^ unique root, and de- 
pends on the k initial values y QI — Note that Y n 

is increasing if any root x^ >1, and also if x^ = 1 
then any terras w^ n 3 ’* x® will increase if ny > 1. There- 
fore, two different solutions to (8) with init.al values y Q , 

z n define e = y„ - z„, and e_ will be bounded for any hX 
o n ■‘n n n 1 

and any initial conditions if and only if all roots of (11) 
satisfy the two conditions above i.e. the stability region 
S = [hX: all roots of (11) are less than one in norm, or on 
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the unit circle and simple]. Figure 1.2 shows the stability 
regions for selected Adams -Bashforth predictors (see Ap- 
pendix A). These are plotted by allowing x = exp(i<j>) for 
some large number of points 0 between 0 and 2 n, and solving 
(11) for hX in each case. 

If the corrector is computed by Q iterations, then (11) 
becomes instead 


y (0) 



k 

2 

j=l 


((a. + hAbj) x k “ j ) 


( 12 ) 




((a* + hXbj ) x k “3) + hXbg y 



s = 1, — ,Q and for each x = exp (i(j»), there will be Q 
values of hX since y^ is a Q degree polynomial in hX. 
However, the lim (Q-*») y^ = y n , the exact solution to 
(11), if (12) converges at all. 

However, problems of interest are not linear homo- 
geneous, and therefore the linear stability analysis gives 
information only about the local behavior of the related 
differential equation 

z' = J n (z-y(t)) + f(y(t), t) 

Z <V = y n 

This behavior could change drastically for even small 
changes in y , and in many cases the eigenvalues of J R are 
not known unless the equations were artificially linearized 
before solving. A more useful stability analysis would try 
to match the stability characteristics of the nonlinear 
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function f(y,t) by proper choice of numerical method which 
most closely approximates those characteristics in important 
areas of the domain of f(y,t). This has been attempted in 
the work reported here. Software has been written which 
will accept a system of up to 20 nonlinear differential 
equations, specified by input equations similar to simula- 
tion languages such as DARE-P [Korn et. al . ] . This will 
output two subroutines. One of these, DIFFUN (T,Y,DY), will 
return in the array DY(*), *=1, — ,m, the derivatives evalu- 
ate at independent variable t=T, dependent variable Y(*). 
The other subroutine, SOL(T, YO, YNEW, INO) will return, in 
YNEW( * ) , a series solution of the equation at t = T > 0 
given the initial conditions y Q in Y0(*), if possible. This 
is described in Chapter 3. 

These routines are then used in the interactive gra- 
phics software described in Chapter 2 which, given M-2 
initial values, searches the remaining 2-dimensional plane 
about an approximate equilibrium point for a connected 
region of initial values having a particular property re- 
lated to stability. These regions can be graphicaJ ly dis- 
played for both the exact and various numerical solutions. 
The routine DIFFUN can also be used by the testing routine 
CLMP described in Chapter 4. This routine, given DIFFUN, 
the initial values, and possibly an inhomogeneous input 
u(t), will display the Fourier amplitude for the first 20 


10 



harmonics for both an "ideal" solution generated by state of 
the art software, and also by a chosen fixed step size 
numerical method. These amplitudes can be compared graphi- 
cally using a bar graph, two solutions can be graphed, and 
the ratios of the Fourier harmoni :s to the input can be 
displayed. In summary, this software can be used to choose 
the most appropriate of available numerical methods by 
comparing stability domains for numerical solutions to the 
same domains for the exact solution. This will insure that 
the numerical solution is stable. Then it is possible to 
verify that decision by testing for accuracy by observing 
appropriate results of a numerical simulation. Each of the 
following chapters outlines one stage in this sequence; the 
first section of each chapter will give some relevant theo- 
retical considerations, the second section will comprise a 
user's manual for that segment of the software, with exam- 
ples. 
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2. STABILITY ANALYSIS OF THE 
NONLINEAR INITIAL VALUE PROBLEM 
The standard linear stability analysis is formula 
specific and describes the behavior of a numerical formula 
applied to the complex test equation 

Y' = 

y<V = V 

Let E n = [ e (t n _ k+1 h . • . / e (t n _ 1 ), e(t n )] be the difference 
sequence for y(t n )-z(t n ) where each solves the test equation 
for a different initial value y Q , z Q . Then e(t n ) = (y Q - 
z 0 )exp(At) is non- increasing in norm for Real (A) < 0. Such 
a condition is called stability. It is desirable for the 
numerical solution y n to be stable if the true solution is 
stable, so for a given stepsize h, one finds all complex A 
such that any numerical sequence E n = [ e n _k+i' • * • ,e n-l' e n^ 
has the property | e^ +1 1 < | j for all i, where 

e i = ^i “ 2 i' t ^ le difference between the numerical solution 
sequences with initial values y Q , z Q . 

For Euler's formula, y n = (1 + hA) y hence the 
numerical solution is stable for 1 1 + hA J £ 1. For multi- 
step formulas, linear stability is characterized by the 
generating polynomials 

k k-i 

r(x) = I a. x K 1 

i=0 1 

k k-i 

s(x) = I b- x K 1 (1) 

i=l 1 
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where r and s have no common divisors. The region of linear 
stability is all hX such that r(x) + hX s(x) has all roots 
inside the unit circle, or on the unit circle and simple. 
For Euler's formula, the r(x) + hX s(x) is (-x+l)+ hX. 
Since this analysis is formula specific, to investigate the 
formula's effect on an actual f(y,t) one considers all the 
eigenvalues X^ of the Jacobian matrix 3f(y,t)/8y; if all hX^ 
are inside the stability region for all y n , t n of interest, 
then the numerical solution will be stable. Insuring that 
such a condition holds is usually not desirable and often 
not possible. 

To develop a stability analysis for nonlinear f(y,t) 
let f(y,t) have the property 

Real < y-z, f(y t)-f(z,t)> < m || y~z|| (2) 

for all t,y, and z of interest. Here <u,v> = y^ Qv for some 
positive definite Hermitian matrix Q, and ||u|| = <u,u>. 

Then for any two solutions y(t), z(t) = y(t)-e(t), e(t) 
satisfies 

de(t)/dt = f(y(t),t) - f(z(t),t) 
and (2) implies that 

d |i e (t)|| 2 /dt = 2 Real < e(t) ,de(t)/dt> < 2\x | e( t)|J 2 and 

thus 

e(t) < exp(pt) e(t Q ) 
which is non-increasing 7 or u < 0 . 

However, (2) is again a condition that cannot be easily 
verified, so a concept relating the true solution sequence 
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Y(t n ) to the computed solution sequence Y n is presented 
here. The following definitions and theorem are helpful. 
They occur in [Dahlquist, 1978] and the theorem proof is 
presented in [Brown, 1979a] and is reproduced in Appendix B. 
Definition - a linear k step formula satisfies 

0 = j! 0 a i + “i f(? n-i' Vjl < 3 > 


Definition - a one leg k step formula corresponding to 
(3) satisfies 


k k 

0 = i a • y ■ + hsf(l/s Z 
j=0 - 1 n 3 j=0 


k 

b • y • , 1/s Z b • t • ) 
3 *n-j’ 7 jt 0 3 n- 3 ' 


(4) 


where s = s(l) and without loss of generality can be set to 
1 by proper scaling of the coefficients. 

Theorem 1 - Let Y n be a sequence which satisfies (4), 
and let Y n = (y n ) be such that 

k 

5 n = £ 0 b j y n -j = s(E) y n (5 > 

where E denotes the back shifting operator. Then Y n satis- 
fies (3). Conversely, if Y n satisfies (3) then there exists 
a sequence Y n such that y n = s (E) y n , and Y n satisfies (4). 

This shows that Y n given by the one-leg k step formula 
will have similar stability properties to its corresponding 
linear k-step sequence Y. In [Dahlquist, 1975] there is 
described a discrete Liapunov function v q ^ ^ which, applied 
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to a sequence Y n , characterizes the stability of that se- 
quence generated by a nonlinear system y' = f(y,t). 

Let 


'e.k.h, <V A <y n+l-i, 


where G is a positive definite, k by k symmetric matrix. The 
structure of G assures that v q ^ ^ iS P os iti ve definite for 



t { 0 }. 


Definition - The (G, k, h) domain of attraction of 
(the numerical solution to) the system is all z Q such that 

* V G,k,h< Z 0> = V G,k,h < Z 1> - V G,k,h < Z 0> <- °' 

where 


Z Q = {z( (l-k)h), . . . ,z(-h),z(0)} , 
Z 1 — { z (( l-k)h Z q , z^} 
in the numerical case and 


Z x = {z( (l-k)h), . . .,z Q ,z(h)} 


for the exact case. 

Definition - The (G,k,h) stability region cf (the num- 
erical solution to) the nonlinear system is all z Q such that 

V G,k,h < Z 0> S inf V G,k,h < Z 0> 

(2 0 c3D) 


where 8D is the boundary of the (G,k,h)-domain of attrac- 
tion. 


This has the following application. Rather than re- 
quiring that Real <y-z, f(t,y)-f(t,z)> < 0, a connected 

subset of initial values y Q is found such that y(h) will be 
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in that subset if y Q was. This is the stability region. 
This insures that the difference y(h)-z(h) is bounded since 
both y(h) and z(h) are in the stability region if y Q and z Q 
were. If f(y,t) is autonomous, y(t n ) will remain in the re- 
gion as n-*». For most well-behaved functions f(y,t), the 
boundary of the region around a stable point can be approxi- 
mated computationally. Once the analytic stability region 
is known, the numerical stability region can be calculated 
using the one-leg k-step method for the same sequence 
[y( (l-k)h), . . . , y(h), y Q ] to get y 1 . The two regions can 
then be compared. 

Analytically, it is possible to form a particular G, 
based on the coefficients of a one-leg k step method, such 
that all numerical sequences based on f(y,t) that satisfy 
(2) will have a stable solution. It was shown in [Liniger 
and Odeh] how to pick G for second order two step formulas, 
second order three-step formulas, and third order three-step 
formulas. 

It is shown below that even an arbitrary choice of the 
positive definite Hermitian matrix G will generate some 
usable results, and theorem 2 demonstrates that using some G 
for a one-leg k step solution y n will generate the same 
stability region for the related solution y n of the linear k 
step formula for a modified 5 . The proof appears in [Brown 
1979a] and in appendix B. 
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Theorem 2 - v q k h = c for t * ie SY 1 * 1111 ®^* 0 posi- 

tive definite matrix G, then there exists a symmetric, posi- 
tive definite matrix 5 , dependent only on G and s(x), such 
that V G k h (? n ) = c, where y n = S(E)y n are the elements of 

and Y_. 
n n 

Sufficient conditions can be developed for the exis- 
tence of the (G,k,h) stability regions based on known tech- 
niques such as Liapunov's direct method [Lenigk] and pro- 
perty (2). An important concept in the development is the 
equilibrium y (t) of the differential function f(y,t). 
While some references define it for an initial value 

y* (t Q ) = 0 

in the space of the dependent variables, this is accom- 
plished by an unnecessary change of variables that could be 
confusing. The important point is that y (t) satisfies 

f(y*(t 0 ),t Q ) be 0 at t Q = 0, 

so that, if f is autonomous', then y(t) is constant, and 

. . . , * 

otherwise, the Taylor series about t Q is given by y (t Q ) + 

t-t Q ) 2 f'(z)/2 and is thus slowly varying and nearly con- 
stant for t near t Q . 

jjf 

Definition - The solution y (t) of y' = f(y,t) for 
y(tg) = y Q , such that f(y Q , t Q ) = 0, is called the equili - 
brium of f(y, t) . 

Definition - The equilibrium of f(y,t) is said to be 
asymptotically stable if there exists a t^ in (t Q , *) such 
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that for every e>0 there exists = d x (e,^) > 0 such that 

if | y 0 ’ yj | < d l' than 

|y(t)-y*(t)|< e 

in [t 1# *]# and there exists a t 2 in (t Q ,») and d 2 (t 2 ) such 
that if |y Q - yj | < d 2 (t 2 ) then 

lim(t+ •) | y(t) - y*(t)| = 0, 
where y(t) satisfies 

Y'(t) = f(y,t), 

y^o* = y<>- 

Theorem 3 [Lenikq] - The equilibrium of f(y,t) is 
asymptotically stable if there exists a function v(y,t) 
which is positive definite in some region D about y*(t Q ) and 
lim v(y,t) = 0 uniformly in t as ||y-y*|| 0 there, and whose 

total derivative is negative definite on D. 

With this background, it can be shown that a well- 
behaved f(y,t) which has a not necessarily unique Liapunov 
function v(y,t) with region D implies the existence of a 
(I,k,h) domain of attraction D' and stability region D" of 
both the exact solution and of a one-leg k step solution 
based on a stable multi-step formula. These two theorems 
are stated and proved in Appendix B and in [Brown, 1979b]. 

In [Dahlquist, 1978] and [Leveque et al] an algorithm 
is presented for calculating a G-stability matrix given any 
A-stable linear multistep method. G is guaranteed to be 
positive definite and symmetric. This code is included in 
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STAN. Some notes on theoretical considerations concerning 
its use are presented here. Since G can only be generated 
for A-stable formulas which include the entire left complex 
half-plane in their stability region, and yet only certain 
implicit multistep methods of order no greater than two can 
be A-stable, the method is not immediately useful for any 
explicit and most implicit formulas. However, by considering 
the modified method 

r (x) = ar(x) + bs(x) 

it 

s (x) = cr(x) + ds(x) 


then 

for 


0 = r*(E)y n + qs* (E)y n 


q = (aA + b)/(cA + d) 
may be A-stable, even though 

0 = r(E)y n + As(E)y n 

is not, for proper choice of a, b, c, and d. 

A further extension looks at 

** * * 
r(x) = r (<t>x) - m(<|> )s (<|>x) 

it it it 

s (x) = S ( <J)X ) 


where 

A A 

m{ <J> ) - min (x = <}>) Real (r (x)/s (x)). 

The coefficients a, b, c, d, and <t> have the following ap- 
plications in the work here. If the multistep formula is 
stable at infinity, then there exists some point m 1 such 
that the linear stability region includes hA such that 
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Real (hA) < Whether or not the formula is stable at in- 
finity, it will contain a region in the left half-plane 
adjacent to 0., and we are interested in the largest stable 
disk, with diameter m 2 , tangent to the imaginary axis at 0.. 
These cases are handled by setting 0-1, and (a, b, c, d) = 
(1,0, 0,1) and finding m^ = m(l) for the stable at infinity 

it it it 

case for r , s . For the disk case, set 0-1, (a,b,c,d) 

= (0, 1, 1, 0 ) , and the diameter is -l/m(0). Figure 2.1 il- 
lustrates this for the 5^ order BDF formula. 

After obtaining G by these techniques, it is guaranteed 
that if Y n -s generated using (4) from 

y' = f(y,t) 

then 

V G,k,h <W S V G,k,h <V 

where 

{ 0 if Mb < m(0) 

0( (l+b(0 ) ( M h-m<0 ) )/(l-b0 ) ( M h-m(0 ) ) ) 1/2 
if m(0) < Mb < m(0) + l/h(0) if Mb > m(0) 

where 

Real < ahf (u(y ) ) + bu(y )-( ahf ( v(y ) ) + bv(y ) ) , u(y )-v(y )> 

iM || u(y)-v(y)|| 2 

and u and v satisfy 

chf(u(y ) ) + d(u(y ) ) = y. 

At A 

The term b(0) depends on s . 
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Figure 2.1. Interpretation of mU) for methods stable 
at Infinity (m. } and not necessarily stabl 
at Infinity (m^). 
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However, we are interested in those initial values 
corresponding to Real < f(y)-f(z), y-z> <0, but computing 
where 

V G,k,h <*i> i V G.k,h<V 

G chosen as above, will show the particular numerical method 
to its best advantage. In the analytic case, G is chosen so 
that 

G = {g kk = 1, g A j = 0 otherwise}. 

Since V(Y) C k h depends on tht oss-product 

<y if yj> s y\ Qyj 

it is helpful to compute Q that is app» . : riate to the two- 
dimensional subspace of the problem being solved. This is 
done by computing the positive definite, symmetric matrix 
that would be used as the Liapunov stability function v(y) 
if only those two variables were involved. This is done by 
decomposing 

f(y,t) = By + g(y, t) 

where B is the numerically differenced Jacobian around t = 

0, y = y* the approximate equilibrium. The Q can be com- 
puted so that 

+ QB = I 

for ^ the transpose operator. 
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STAN USER MANUAL 


An interactive graphics program called STability ANaly- 
sis (STAN) is stored in object form on disk and can be 
linked with two subroutines SOL (T, YO, YNEW, IND) and DIFFUN 
(T, Y,DY) . These subroutines can be written by the user or 
be produced by the PASCAL program SERIES described in Chap- 
ter 3. DIFFUN returns in the N-dimensioned array DY the 
derivative f(y,t) given t = T, Y=y(t) an array of depen- 
dent variable values. SOL returns in the N-dimensional 
array YNEW the solution to y'=f(y,t) at T=t, given YO = y ( 0 ) 
an N-dimensioned array of initial values at t=0. SERIES 
uses recursive power series techniques for this and sets IND 
> 0 if the radius of convergence is not (-1,1); if the ana- 
lytic solution is known it can be programmed by the user. 
STAN can also work with only a dummy routine SOL (T,YO, 

YNEW, IND) which assigns YNEW(I) = YO(I), 1=1,..., N, if the 
analytic stability properties are never requested. 

The program uses prompts to guide the user through its 
execution. When the possible commands, usually strings of 
no more than 10 characters, are not obvious, the command 
HELP will list them. Only enough of the command to estab- 
lish its unique identity need be given. Many of the com- 
mands are not needed for simple jobs because the execution 
of the GO command will automatically prompt the user to 
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provide the necessary information. The commands are avail- 
able to allow the user to reset his system for more compli- 
cated jobs. Default values are provided for almost every 
program parameter. A typical load under the NOS operating 
system, given the object code of STAN and the subroutines 
SOL and DIFFUN from SERIES, follows - 
ATTACH, STAN. 

ATTACH, SOL. 

ATTACH, DIFFUN. 

FTN , I =SOL , B=X . 

FTN, I=DIFFUN,B=Y. 

LIBRARY( IMSLLIB, PLOTIO ) 

LOAD ( STAN , X , Y ) 

NOGO , Z . 

PLOTIO and IMSLLIB are the TEKTRONIX graphics library and 
the IMSL library, respectively. 2 is the resulting load 
module, to be executed. 

The available commands are of three types - those that 
describe the system of equations being investigated; those 
describing the solution method being used, including the 
analytic solution; and those describing the desired display. 
These will be described here. The first group includes 
INIT, CENTER, SYSTEM, NEQ, and NORM. A two-dimensional 
subsystem of N-differential equations will be investigated 
by holding N-2 variables constant at t~0, and varying the 
other 2 "active" variables about an approximate equilibrium 
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point CENTER where the derivatives of the two active vari- 
ables are nearly zero. Such a point must be known in ad- 
vance from analytic considerations, but the equilibrium can 
be improved by using Mueller's method [Conte, de Boor] for 
finding the root of 2-dimensional nonlinear function. The 
commands are - 

NEQ - specify N, the dimension of the system. Default is 
two. 

SYSTEM - specify il, i2 the two dependent variables to be 
investigated. Default is il = 1, i2 = 2. 

INIT - set the values the N-dependent variables take at t = 

0 . 

Default is 0. 

CENTER - set the equilibrium point CENTER of the two- 
dimensional subspace where f(y,t) = 0. When the program 
prompts IMPROVE CENTER, . . . , a reply of YES will use the 
initial value given and try to solve 

f il ( y 't ) =0 

and 

f i2 (y,t) = 0. 

NORM - find Q such that Q+QB=I where B is the 2 by 2 
Jacobian of the system with respect to y^, y^ 2 - Whenever 
SYSTEM is called after the initial GO command, this is done 
automatically. Calling NORM is usually not necessary unless 
the stepsizc H has been changed, or the initial values are 
drastically changed. B is computed by differencing. 
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Commands describing the solution being used are NSTEPS, 
STEPS I ZE, and TYPE. These commands tell whether the ana- 
lytic solution or a k-step numerical method is being used, 
and what stepsize H is employed. Since the desired output 
is either D' or D", initial values of y^ and y^ 2 such that 
the function V(Y) G ^ ^ is decreasing, these commands are 
used to set k and h. Note that one should compare the 
analytic D' and D 1 ' for the same k as the numerical method, 
since D’ and D M are likely to be smaller as KH grows larger 
KH < 1 is required if SOL is generated by a power series 
expansion. G is picked automatically for the various numer- 
ical methods, G = [g kk = 1, g^j = 0 otherwise] for the exact 
case, so the interior of D' satisfies 

<y(kh) ,y(kh)> < < y( (k-1 )h) ,y( (k-1 )h)> 
for the exact solution y(t). 

The commands are - 

NSTEPS - set K, the number of steps in the numerical methods 
being compared. When K is changed, D' and D' 1 for the exact 
case should be recomputed. Default is 1, maximum is 10. 
STEPSIZE - H, the stepsize in the independent variable, is 
set. This is normally set by GO the first time, and need 
not be changed unless D' is needed for various stepsizes. 
TYPE - EXACT or NUMERICAL, sets indicators so that when the 
GO command is given a K step numerical method and correspon- 
ding matrix G are chosen. 
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The commands that describe the display are NPOINTS, 
SYMMETRIC, PARAM, and REPEAT. 

NPOINTS - sets NP, the NP rays spaced at equal angular dis- 
placement about CENTER in R 2 are generated to find the 
boundaries of D' and D M . 3 < NP < 80, default 13. 

SYMMETR T C - if the system is known to be symmetric, then the 
NP plotted points are concentrated in only half of R 2 . When 
the system prompts is answered by 1, this means the system 
is symmetric about 0 (the upper half-plane is graphed); 
2 means symmetry about y^ = 0 (the right half plane is 
graphed); 0 cancels the symmetric display. 

PARAM - display most of the parameters that have been set. 
REPEAT - repeat the last graph, allowing display limits to 
be set. 

The command GO initiates calculation of D' and D’ ' . On 
the first GO call, the user must supply H, then Q is com- 
puted and listed. If this is a numerical test, one of the 
numerical methods in Appendix A must be chosen, or else a k 
step predictor corrector (PC) is input by the user. When 
asked for the predictor, enter 

a 0 , a^, • • • * a^, bQ— 0, b^, . . « , bj^. 

If the number of corrector iterations is given as 0, no 
corrector need be entered. Otherwise, give 



The program computes G, depending on whether the corrector 
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is known to be stable at infinity or not. If the approx- 
imate equilibrium is not inside D, the Liapunov stability 
region, then the message INVERSE REG. ON will be printed; try 
a new CENTER or a different subsystem. If the system is too 
irregular or too nonlinear, then the message UNABLE TO... 
and a PAUSE will be printed. A carriage return will cause 
the program to return to the beginning command sequence. 
Choice of a different subsystem or rewriting the system of 
equations may solve this problem. 

After D', the domain of attraction, is computed, it can 
be displayed. On the first GO call, the left, right, bot- 
tom, and top values in D' are printed, and the user can 
choose his display coordinates. Thereafter, these coordi- 
nates can be changed or left alone. After choosing whether 
to display D', one can compute and display D ,! in a similar 
fashion with the default coordinates being those chosen for 
D’. In many problems, D* ' may be so close to D* that only 
D' need be displayed. This w>'ll save a great deal of compu- 
tation. After this, the user is returned to the beginning 
command sequence, but now all parameters have been initia- 
lized and, after changing any of them, a new GO will proceed 
faster than the first. 
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EXAMPLE 


The following set of equations form a subsystem that 
models the longitudinal stability of the F4 aircraft [Stein- 
metz et al] - 

A = ATAN ( W/U ) 

V 2 = U 2 + W 2 

U' = -C1*SIN 0 - W*Q + (Cl + C2*A)V 2 + C3*Q*V 
W* = Cl*COS 0 + U*Q + ( C4+C5 *D+C6 * A ) V 2 + C7*Q*V 
Q» = (C8+C9*(A+D))V 2 + C10*W'/U + C11*Q*V 

0’ = Q 

D = 2*T* . 174533 if 0 < T < 0.5 
0.1734533 if 0.5 < T 

where U is the horizontal velocity along the aircraft body, 
W is the vertical velocity perpendicular to U, A is angle of 
attack, 0 is pitch, Q is rate of change of pitch, and D is 
the driving function, the stabilator deflection angle. These 
equations will be put in a usable form for input to SERIES 
and a sample run of STAN will illustrate t} e various fea- 
tures. The constants Ci appear in the program listing in 
Figure 2.2. 
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Figure 2.2a Input to Series 


/* C=32.16, Cl=- . 5Q18521E-9 , C2=. 1192125E-5, C3=- .42110675 
E-3, C4=- . 14598254E-3 , C5=.46345531E-5, 
C6=-.45006065E-5, C7=-.6971899E-3, C8=-.587215E-8, 
C9--.73872E-6, C10=. 3800645E-5, Cll=- . 7422436 E-2*/ 

U . =-C*S IN ( THETA ) -W*Q+ ( C1+C2 *A ) * ( U* *2 +W* *2 ) 

+C3*Q*(U**2+W**2 )**.5; 


W . =C*COS ( THETA ) +U*Q+ ( C4+C5*D+C6*A ) * (U**2+W**2 ) 
+C7*Q*(U**2+W**2)**.5; 


Q.=(U**2+W**2)*(C8+C9*(A+D) )+(C10*( (C*COS(THETA) 
+U*Q+ ( C4+C5 *D+C6 *A ) * ( U* *2+W* *2 ) 
+C7*Q*(U**2+W**2 ) ** . 5 )/U)+Cll*Q ) 

*(U**2+W**2 )** .5; 

THETA. =Q; 


A. =(U*( C*COS ( THETA ) +U*Q+ ( C4+C5 *D+C6 *A ) * ( U* *2+W* *2 ) 
+C7*Q*(U**2+W**2 )** . 5 )-W*( -C*SIN(THETA) 
-W*Q(C1+C2*A)*(U**2+W**2) 

+C3*Q*(U**2+W**2 )**. 5 ) )/(U**2+W**2 ) ; 

D.=. 349066; ; 
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Figure 2.2b. 

DIFFUN, the output from SERIES used by STAN and CLMP. 

SUBROUTINE DIFFUN(T, Y,DY) 

DIMENSION DY(20), Y(20) 

DATA C/32.16/, C1/-.5018521E-9/, C2/. 1192125E-5/, 

C3/-.421106759E-3/, +C8/-.587215E-8/, C9/-.73872E-6/, 

C10/.3800645E-5/, C11/-.7422436E-2/ 

DY(1 )=-C*SIN(Y(4) )-Y(2 )*Y(3 )+(Cl+C2*Y(5 ) )*(Y(1 )**2+Y(2 )**2 ) 
+C3*Y(+3 )*(Y( 1 )**2+Y(2 )**2 )** . 5 

DY(2 )=C*C0S(Y(4) )+Y( 1 )*Y(3 )+C4+C5*Y(6)+C6*Y(5 ) )*(Y(1 )**2 
+Y(2 )**2+)+C7*Y(3 )*(Y(1)**2+Y(2 )**2)** .5 

DY( 3 )=(Y( 1 )**2+Y(2 )**2 )*(C8+C9*(Y(5)+Y(6 ) ) ) 

+(C10*( (C*COS(Y(4) )+Y(+l)*Y(3) 
+(C4+C5*)(6)+C6*Y(5))*(Y(1)**2+Y(2)**2) 

+C7*Y(3 )*(Y(1 )**2+Y(+2 )**2 )** . 5 )/Y(l ) 

+C11*Y(3 ) )*(Y(1)**2+Y(2 )**2 )**.5 

DY ( 4 ) -Y ( 3 ) 

DY(5)=(Y(l)*(C*COS(Y(4))+Y(l)*Y(3)+(C4+C5*)(6) 

+C6*6(5 ) )*(Y(1)**2+Y(2)**2 )+C&*Y(3 )*(Y(2 )**2+Y 
(+2 )**2 )** . 5 )/Y(2 ) )+Cll*Y(3 ) )*(Y(2 )**2+Y(2 )**2 )** . 5 

DY(4)=Y(3 ) 

DY(5) + (Y(1)*( C*COS (Y(4))+Y(l)*Y(3)+( C4+C5 *Y ( 6 ) +C6 *Y( 5 ) ) 

*(Y( 1)**2+Y(2 )**2(+C7*Y(3 )*(Y(1 )**2+Y(2 )**2 ) 
**.5)Y(2)*(-C*SIN(Y(4))-Y)2*Y(+3) 

+(C1+C2*Y(5) )*(Y(1)**2+Y(2)**2) 
+C3*Y(3)*(5f(l)**2+Y(2)**2)**.5) )/ 

+(Y( 1)**2+Y(2 )**2 ) 

DY(6)=. 349066 

RETURN 

END 
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Figure 2.2C - SOL, used by STAN. 

SUBROUTINE SOL(T, YO, YNEW, IND) 

DIMENSION Y0(20) , YNEW(20), 2Z2B(20), TBC(20), DW (20), 

TB7(20), TBD(20), +D(20), TB8(20), TBE(20), TC7(20), 

TDA(20 ) , TBF(20 ) , TCD(20), TC8(20) ,TBG+(20) , TCE(20), 

TC9(20), TBH(20 ) , TBI (20 ) , TBJ(20), TDF(20), TBK(20), 

TCI+(20) , TDG(20 ) , TBL(20), TCJ(20), DA(20), TBM(20), 
TCK(20 ) , TBN( 20 ) , TCL+20, DTHETA(20) , DD(20), Q(20), 

TBR(20 ) , TBS (20 ) , TBT(20), TEN(20), U+(20), TBU(20), 

TEO(20) , TBV(20) , W(20), TBW(20), TBX(20), THETA(20) , 
TBY+20), TES(20 ) , TBZ(20), TBO(20), TB1(20), DQ(20), 
TB2(20), DU(20 ) , TBA+(20), A(20), TBB(20), TB6(20), 

TD2(20 ) 

DATA C/32.16/, C1/-.5018521E-9/, C2/. 1192125E-5/, 
C3/-.42110675E-3/, C+4/- • 14598254E-3/, 

C5/.46345531E-5/, C6/- .45006065E-5/, 

C7/- . 6971899+E-3/, C8/- . 587215E-8/, C9/- . 73872E-6/, 

CIO/. 3800645E-5/, Cll/- . 74224+36E-2/ 

EPS=1 . 0E-10 
U(l)=YO(l) 

W(1 )=YO(2 ) 

Q( 1 )=YO(3 ) 

THETA(l)=YO(4) 

A ( 1 ) =YO ( 5 ) =ATAN2 ( YO ( 2 ) , YO(l)) 

D(l)=YO(6) 

IND=0 

DO 1 111=1,19 
NIII=III 

1111=111-1 

IF(III.EQ.l) GO TO 100 
TBB( III )=0. TBA( III )+0. 

DO 101 JJJ=1,IIU 

TBA( III )=TBA( III )+TBB( JJJ) *( I II-JJJ)*THETA( I II-JJJ+1 ) 

101 TBB( III )+TBB( III )-TBB( III )-TBA( JJJ)*( III-JJJ) *THETA 
(III-JJJ+1) 

TBA( III )=TBA( 1 1 1 )/ ( 1 1 1 — 1 . ) 

TBB( III )=TBB( III )/( III-l . ) 

GO TO 102 

100 TBA (III) =S IN ( THETA (III)) 

TBB (III) +COS ( THETA (III)) 

102 CONTINUE 

TBC (III) =TBA ( 1 1 1 ) *C 
TBD (III ) =-TBC ) 1 1 1 ) 

TBD( III) =0 . Do 103 JJJ=1,III 
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103 


Do 103 JJJ=1,III 

TBE( III )=TBE(III )+W( JJJ) *Q( III-JJJ+1) 

TBF( III )=TBD( III )-TBE( III ) 

TBG(III )=A(III )*C2 
TBH( III )=TBF( III ) 

IF ( I I I . EQ. 1 )TBH( I I I )=C1+TBG( III) 

TBI (III ) s O. 

Do 104 JJJ=1,III 

104 TBI ( III )=TBI( III )+U(JJJ)*U( III-JJJ+1) 

TB J ( I I I ) =0 . 

Do 105 JJJ=1,III 

105 TB J (III) =TB J ( I I I ) +W ( JJJ )*W( III-JJJ+1) 

TBK( III )=TBI (III )+TBJ (III )TBL( III )=0. 

DO 106 JJJ=1, III 

106 TBL( III )=TBL( I II ) +TBH ( J J J ) *TBK ( III-JJJ+1) 

TBM( III )=TBF( III )+TBL( III ) 

TBN(III)=Q(III)*C3 

IF(III.EQ.l) Go to 107 
TBR( III ) = 0. 

Do 108 JJJ=1,III1 

108 TBR( III )=TBR( III )+(((. 5+1. )*( III-JJJ) ) 

/( III-l . )-l. )*TBR(JJJ)*TBK+( III-JJJ+1) 

TBR(III )=TBR(III )/TBK(l) 

Go to 110 

107 TBR(III )=TBK(III )**.5 

110 CONTINUE 
TBS ( I I I ) =0 . 

Do 111 JJJ=1,III 

111 TBS ( I I I ) =TBS (III) +TBN ( J J J ) *TBR ( I I I - J J J+ 1 ) 

TBT (III) =TBM (III) +TBS (III) 

DU( III )=TBT( III ) 

U( I II+l )=DU( III )/FLOAT( III) 

IF (III.EQ.l) GO to 112 
TBU ( III )=0. 

TBV( III )=0. 

Do 113 JJJ=l t IIIl. 

TBV( III )=TBV( III )+TBU(JJJ)*( III-JJJ )*THETA( III-JJJ+1) 

113 TBU ( I I I )=TBU( III) -TBV( JJJ ) * ( I I I -JJJ ) * THETA ( I I I-JJJ+1 ) 
TBV( III )=TBV( I I I )/( II 1-1 . ) 

TBU(III )=TBU(III )/( III— 1. ) 

Go to 114 

112 TBV( III )=SIN(THETA( III)) 

TBU (III) =COS ( THETA( III)} 

114 CONTINUE 

TBW (III) =TBU (III) *C 
TBX( II I )=0. 

DO 115 JJJ=1,III. 

115 TBX( III )=TBX( III) +U (JJJ) *Q (III-JJJ + 1) 

TBY ( 1 1 1 )=TBW( 1 1 1 )+TBX( III) 

TBZ( III) = D (III )*C5 
TB0 (III) =TBZ (III) 
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I 

| If (III.EQ.l )TBO( III ) S C4+TBZ( III) 

! TB1( III ) S A( III )*C6 

TB2 (III )=TBO( III )+TBl( III ) 

, TB6(III ) s 0. 

Do 116 JJJ»1,III. 

116 TB6( III ) S TB6( III )+TB2( JJJ) *TBK( I II-JJJ+1) 
TB7 (III) =TBY (III) +TB6 (III) 

l TB8( III )=Q( III )*C7 

i, TCD( III )=0. 

117 TCD ( I I I ) =TCD( III) +TB8 ( JJJ ) *TBR ( I II-JJJ+1 ) 

t TCE( III )=TB7( III )+TCD( III ) 

DW( III )=TCE( III ) 

W( III+1)=DW( III )/FLOAT( III ) 
TCI(III)=A(III)+D(III) 

TCJ( III ) S TCI( III )*C9 

* TCK( III )=TCJ( III ) 

IF(III.EQ.l )TCK( III )=C8+TCJ( III) 

\ TCL( III )=0. 

I DO 118 JJJ=1,III 

118 TCL( III ) =TCL( III) +TBK( JJJ ) *TCK( I I I - JJJ+1 ) 

, IF (III.EQ.l) GO TO 120 

TC7(III)=TCE(III)-TC7(1)*U(III) 

* I F ( I I I . EQ . 2 ) GO TO 121 

DO 122 JJJ=2, III1 

{ 122 TC7 (III )=TC7 (III )-TC7( JJJ) *U( II I- JJJ+1 ) 

L 121 TC7( III )=TC7( III )/U( 1) 

GO TO 123 

i 120 TC7 < 1 1 1 ) +TCE (III) /U (III) 

123 CONTINUE 

TC8( III )=TC7( III )*C10 
TC9(III )=Q(III )*C11 

1 TDA(III)=TC8(III)+TC9(III) 

■ TDF( III )=0. 

DO 124 JJJ=1,III 

I' 124 TDF( III )=TDF( III )+TDA( JJJ)*TBR( I II-JJJ+1 ) 

! TDG(III )=TCL(III )+TDF(III ) 

DQ( III )=TDG( III) 

j Q( imi )=DQ( III )/FLOAT( II I ) 

! DTHETA( III ) S Q( III ) 

THETA ( I I 1+1 )=DTHETA( I I I ) /FLOAT ( III) 
TD2(III )=0. 

! DO 125 JJJ=1,III 

! 125 TD2 (III )=TD2 (III )+U( JJJ)*TCE{ III -JJJ+1) 

TEN ( I I I ) =0 . 

r DO 126 JJJ=1,HI 

126 TEN( III )=TEN( III )+W( JJJ )*TBT( III JJJ+1) 
TEO (III) =TD2 (III) -TEN (III) 

IF (III.EQ.l) GO TO 127 

TES( III )=TEO( III )-TES( 1 )*TBK( III) 

IF( III .EQ.2 ) GO TO 128 
DO 130 JJJ=2,III1 
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130 TES( III )=TES( III )-TES( JJJ)*TBK( III-JJJ+1 ) 
128 TES( III )=TES( III )/TBK(l) 

GO TO 131 

127 TES(III )=TEO(III )/TBK(III ) 

131 CONTINUE 

DA( III ) s TES( III ) 

A( I II+l )=DA( III )/FLOAT( III ) 

IF( III .GT. 1 )DD( III) —0 . 

DD(1)=. 349066 

D( III+l )=DD( III )/FLOAT( III ) 

IF ( III .LT.4)GO TO 1 
1111=111+1 
ZZ221=0 . 

22222=0. 

DO 132 JJJ=1,III1 
ZZZZ1 S ZZZZ1+U ( JJ J ) 

IF( JJJ.LT. III-4) GO TO 132 
ZZZZ2=ZZZZ2+ABS(U{JJJ) ) 

132 CONTINUE 
ZZZZ1=EPS*(ABS(ZZZZ1)+1. ) 

I F ( ZZZZ2 . GT . ZZZZ1 ) GO TO 1 
ZZZZ1=0. 

ZZZ22=0 . 

DO 133 JJJ=1,III1 
ZZZZ1=ZZ2Z1+W( JJJ ) 

IF( JJJ.LT. Ill —4 ) GO TO 133 
ZZZZ2=ZZZZ2+ABS(W( JJJ) ) 

133 CONTINUE 
ZZZZ1=EPS*(ABS(ZZZZ1 )+l . ) 

I F ( ZZZZ2 . GT . ZZZZ1 ) GO TO 1 
ZZZZ1=0 . 

ZZZZ2=0. 

DO 134 JJJ=1,III1 
ZZZZ1=ZZZZ1+Q( JJJ ) 

I F( JJJ.LT. I I 1-4) GO TO 134 
ZZZZ2=ZZZZ2+ABS(Q(JJJ) ) 

134 CONTINUE 
ZZZZ1=EPS*(ABS(ZZZZ1)+1. ) 

I F ( ZZZZ2 . GT . ZZZZ1 ) GO TO 1 
ZZZZ1=0. 

ZZZZ2=0. 

DO 135 JJJ=1,III1 
ZZZZ1=ZZZZ1+THETA( JJ J ) 

IF (JJJ.LT. III-4) GO TO 135 
ZZZZ2=ZZZZ2+ABS( THETA ( JJJ) ) 

135 CONTINUE 
ZZZZ1=EPS*(ABS(ZZZZ1 )+l . ) 

I F ( Z2ZZ2 . GT . ZZZZ1 ) GO TO 1 
ZZZZ1=0. 

ZZZZ2=C . 
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DO 136 JJJ=1,III1 
ZZZZ1=ZZZZ1+A(JJJ) 

IF ( JJJ.LT. III-4) GO TO 136 
ZZZZ2=ZZZZ2+ABS ( A( JJJ ) ) 

136 CONTINUE 

ZZZZ1=EPS*(ABS(ZZZZ1)+1. ) 

IF ( Z2ZZ2 . GT . ZZZZ1 ) GO TO 1 
ZZZZ1=0. 

ZZZZ2=0 . 

DO 137 JJJ+1,III1 
ZZZZ1=ZZZZ1+D( JJJ) 

IF( JJJ.LT. III-4) GO TO 137 
ZZZZ2=ZZZZ2+ABS(D(JJJ) ) 

137 CONTINUE 
ZZZZ1=EFS*(ABS(ZZZZ1 )+l. ) 

I F ( ZZZZ2 . GT . ZZZZ1 ) GO TO 1 
GO TO 2 

1 CONTINUE 

2 CONTINUE 

DO 138 JJJ=1,NIII 

IF(ABS(TBK( JJJ) ) . LT.EPS) GO TO 138 

KKK-JJJ 

GO TO 140 

138 CONTINUE 

140 ZZZZ1-0. 

KKK1=KKK+1 

DO 141 JJJ=KKK1,NIII 

141 ZZZZ1=ZZZZ1+ABS(TBK(JJJ)) 
IF(ZZZZ1/ABS(TBK(KKK) ) .GE. 1 ) IND=IND+1 
DO 142 JJJ=1,NIII 

IF (ABS (U( JJJ)). LT.EPS) GO TO 142 

KKK=JJJ 

GO TO 143 

142 CONTINUE 

143 ZZZZ1=0. 

KKK1=KKK+1 

DO 144 JJJ=KKK1,NIII 

144 ZZZZ1=ZZZZ1+ABS(U(JJJ)) 
IF(ZZZZ1/ABS(U(KKK) ) .GE. 1 ) IND S IND+1 
DO 145 JJJ=1,NIII 

IF(ABS(TBK( JJJ)). LT.EPS) GO TO 145 

KKK=JJJ 

GO TO 146 

145 CONTINUE 146 
ZZZZ1=0 . 

KKK1=KKK+1 

DO 147 JJJ=KKK1,NIII 
147 ZZZZ1=ZZZZ1+ABS(TBK(JJJ)) 

I F ( ZZZZ 1 /ABS ( TBK ( KKK ) ) . GE . 1 ) I ND= I ND+1 
NIII+NIII+1 
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ZZZB(1 )=U(NIII ) 

DO 148 JJJ=2,NIII 

148 ZZZB( JJJ )=U(NI I I - JJJ+1 ) +T*ZZZB( JJJ-1 ) 
YNEW ( 1 ) =ZZZB ( NI 1 1 ) 

ZZZB(l)=W(NIIi) 

DO 150 JJJ=2,NIII 

150 ZZZB ( J J J ) =W(NI I I -JJJ+1 ) + ( T*ZZZB ( JJJ-1 ) 
YNEW(2)=ZZZB(NIII) 

ZZZB91 )=Q(NIII ) 

DO 151 JJJ=2,NIII 

151 ZZZB ( J J J ) =Q ( NI I I -JJJ+1 ) +T*ZZZB ( JJJ-1 ) 
YNEW(3)+ZZZB(NIII) 

ZZZB( 1 )=THETA(NI II) 

DO 152 JJJ=2, Nil I 

152 ZZZ % J J J ) =THETA ( NI I I -JJJ+1 ) +T*ZZZB (JJJ-1 ) 
YNEW(4)=ZZZB(NIII) 

ZZZB(1)=A(NIII) 

DO 153 JJJ=2,NIII 

153 ZZZB ( J J J ) =A ( NI I I -JJJ+1 ) +T*ZZZB ( JJJ-1 ) 
YNEW( 5 )=ZZZB(NIII ) 

ZZZB(1)=D(NIII ) 

DO 154 JJJ=2,NIII 

154 ZZZB ( J J J ) =D (NI I I -JJJ+1 ) +T*ZZZB ( JJJ-1 ) 
YNEW(6 )=ZZZB(NIII ) 

RETURN 



Since SERIES only accepts input with derivatives on the left 
of the =, and a limited subset of FORTRAN functions of 
derivative free variables on the right, the equations for A, 

V 2 , V, and D must be changed. For example, since A = ATAN(W/U), 
the differential equation A' = (UW' - WU')/(U 2 ) could be 
added, but then the entire expression for W and U' would 
have to replace these values (the actual output SOL of 
SERIES will have these common subexpressions eliminated; the 
restriction was imposed to remove the problem of sorting the 
input to SERIES. Sorting often is the cause of errors in 
simulation languages such as ACSL and DARE.) Further, the 
original formulation includes the approximation 

A' = W’/U 

so this slightly simpler expression can be used to change 
from an algebraic to a differential equation. Since A(0) is 
a function of W(0) and U(0) and SERIES cannot handle alge- 
braic constraints, the output of SERIES in Figure 2.2 has 
been modified by the lines Y(5) = ATAN2(Y(2)/Y(1) ) in DIFFUN 
and A ( 1 ) =YO ( 5 ) =ATAN2 ( YO ( 2 )/Y0 ( 1 ) ) in SOL to provide this 
algebraic constraint. 

Since the stepsize H usually used in a real time simu- 
lation package is h = 0.032 and we are limited to k=10 
steps, stabilator deflections occurring after 0.32 sec will 
never occur, so the proper equation for D is 

D' = 2 .*. 174533, 

D( 0 ) = 0 
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V 2 and V are replaced by (U*U+W*W) and SQRT(U*U+W*W) where 
they occur. The resulting system of equations becomes 
U* = -C1*SIN o - W*Q+(C1+C2*A) (U*U+W*W) + C3*Q*SQRT(U*U+W*W) 
W' = C1*C0S 0 + U*Q+(C4+C5*D+C6*A) (U*U+W*W) 

+ C7*Q*SQRT(U*U+W*W) 

Q' = (U*U+W*W) (C8+C9*(A+D) ) + C10*(C1*COS 4> + 
/U*Q+(C4+C5*D+C6*A) (U*U+W*W) )/U+ll*Q*SQRT(U*U+W*W)/U 

<t>' = Q 

D' = 0.349066 

A' = (Cl*COS 0 + (C4+C5*D+C6*A) (U*U+W*W) + 
C7*Q*SQRT(U*U+W*W) )/U 

Initial conditions at t=0 are S = 0, A = ATAN(W,U), = 
5.3° = 0.093 radians, Q = 0. Typical initial conditions for 
U and W are U = 660.18167, W = 5.74626. 

Fig 2.3 shows an extensive terminal session to deter- 
mine initial conditions and best numerical method to be used 
with these equations for explicit multistep or Runge-Kutta 

methods for real-time simulation. Sequence A shows initial- 

. . . . ... * 

ization, choice of approximate equilibrium y = (271,360) 

from an initial guess of U=660. 18167 and W=5. 74626. This 
corresponds to the plane going slow and climbing, and is 
inside the Liapunov stability region. Computation of the 
Liapunov norm matrix shows that only b^ is not almost zero, 
so this means tha W is not sensitive to changes in U or W, 
and that U is not sensitive to changes in W, either. Thus, 
<y,y> = U 2 . Both D' and D M are disp^yed, but since D' ' is 


r 
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proportional to D' only D' will be computed in the fol- 
lowing. 

At B, numerical analysis is chosen, and D' for the AB1 
(Euler) method cannot be plotted. Neither of the two Runge- 
Kutta methods work either, since their stability regions are 
too small, and incidentally identical. 

At C, the number of steps is increased to 2, and the 
exact region D' is computed for k=2. As expected, this 
region is somewhat smaller, at least in the W coordinate. 
At D, the 2 -step Adams predictor is analyzed, and at E, the 
shifted trapezoidal rule. Interestingly, this last has a 
much larger D' than either AB2 or the exact D', but this is 
not an advantage since we want a close match, not the lar- 
gest region. 

At F, a study of the analytic D' for H = 1/24, and 
0.016, shows that, as expected, the larger H is, the smaller 
D' is. 

We can now conclude that a two-step method is needed, 
and at G we look at other subsystems to pick which one. The 
subsystem 2,3 again is only sensitive to its second variable 
Q, but the 1,3 subsystem shows a definite dependence between 
U, the horizontal velocity, and Q, rate of change of pitch. 

At H, we investigate AB2 and shifted trapezoid, and 
conclude that AB2 most closely follows the exact case, since 
shifted trapezoid slows down even when the nose is dropping 
fast ( Q< 0 ) to a greater extent than is the real case. 
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DISPLAY DOMAIN OF ATTRACTION 

? YES 

COMPUTED LIMITS ARE: 

- 23688.9602868 47506 . - 29258.46196759 27000.33534529 

CHANGE DISPLAY LIMITS 



*• *’ 
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COMPUTED LIMITS ARE: 

25.671 19231 207 15496.18167 - 184.4512589434 148.5030431346 

CHANGE DISPLAY LIMITS 




3. NUMERICAL SOLUTIONS 


The analytic solution to the problem 

y 1 “ f<y,t) (1) 

y(0) = y 0 

is attacked numerically for a restricted set of input func- 
tions f using power series methods. That is we represent 
each of the y(i), i=l, n for n <= 20 by a power series 
expanded about t=0, and given t and y Q the value of y(i) can 
be calculated. 

The general method of solution is to use recurrence 
relations to calculate successive terms of the power series. 
Algorithms for addition, subtraction, and multiplication of 
power series are well known. Recurrence relations for 
division and exponentiation can be derived which calculate 

4.L 

the n tn term of the result using only the first n terms of 
the operand or operands and the first n-1 terms of the 
result [Knuth] . The same type of recurrence relations can 
also be derived for sin, cosine, natural logarithhm, and 
exponential of a power series [Gibbons] along with many 
other functions. The computation of the first n terms of 
any of the operations or functions mentioned above can be 
accomplished in 0(n 2 ) or less multiplications and divisions. 
In addition to this restricted set of inputs, any function 
which can be defined as an initial value problem with ini- 
tial conditions at t=0 using the given set of functions and 



operations can be added to the set of equations and conse- 
quently used as part of the input. 

As an example of a recurrence relation, let 

A = la, t 1 ’ 1 , B = Z b, t 1-1 
i=0 1 i=o 

and consider 

B = EXP (A) (2) 

differentiation of both sides with respect to t yields 

B' = EXP (A) * A' 
using (2) 

B* = B * A' 

OB Ob 00 • 

Z nb_ t n ’ x = l b t n * X * Z i a. t 1 * 1 

n=l n n=l n n*l 1 

Z nb t n-1 = Z (Z 11 ’ 1 b t (n-i) a t ) t n ’ x 
n=l n n=l i=0 1 n 1 

equating the coefficients of t n “^ 

nb n = na n b o + + a l b n-l' n - 1 ( 3 ) 

b Q = EXP (a Q ) 

giving us our recurrence relation. 

The general form of the input is: 

/* constantl = valuel, constantn = valuen */ 
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variable!.. = equationl; 

i (4) 

* 

variablem. = equationm; ; 

Constantl, . .., constantn are unique variable names and 
valuel , . . . , valuen are numbers . Derivatives with respect 

to time, represented by variable., can appear explicitly 
only on the left-hand side of the equation. If derivatives 
were allowed to appear on the right-hand side the equations 
would have to be ordered; however, by allowing derivatives 
only on the left-hand side no ordering is necessary by 
either the program or user. 

The input is read as a character string from the data 
file EQIN and assembled into lexical items using a scanner. 
These lexical items are then sent to a recursive descent 
parser which enters them into a symbol table, and generates 
quadruples as the output of each operation, i.e. a*b would 
become *a b tl where the result of a*b would be assigned to 
tl. Finally the codetable is examined to eliminate common 
subexpressions. It is from this optimized code table that 
the recurrence relations are generated, [Gries], [Barton et 
al]. 

The output from the program consists of two FORTRAN 
subroutines. The first is subroutine DIFFUN which is used 
in calculation of the numerical solution by other parts of 
this package and is described in more detail in Chapters 2 
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and 4. The second output is subroutine SOL, the subroutine 
which attempts to solve the system of equations to stated 
accuracy. The constant definition is transformed into a data 
statement and the quadruples generated by the parser are 
used to generate the recurrence relations. Enough terms (up 
to 20) are calculated by SOL so that all the input variables 
satisfy: 

n n 

I ABS(Y(i ) ) < EPS*(ABS(I Y(i) + 1.)) 
i=n-4 " i=0 

as used in [Gibbons). Jn general this condition states that 
the last five terms are negligible compared to the sum of 
the series. Obviously any alternating series will be con- 
vergent using this criteria and for an arbitrary power 
series this test provides some assurance that the rebt of 
the terms of the series can be safely neglected since values 
of t only between -1 and 1 are being considered. 



SERIES USER MANUAL 


The purpose of this part of the software package is to 
generate the FORTRAN subroutines necessary to solve the 
input system of first order ordinary differential equations 
both analytically and numerically. As stated in (4) the 
general form of the input is: 

/*constantl = valuel, . .., constantn - valuen */ 
variablel. = equationl; 

variablem. = equationm;;. 

The first part of the input is the constant definition 
section. Any number of constants can be defined, and these 
will appear in data statements in both subroutines SOL and 
DIFFUN. The second part of the input is the specification 
of the differential equations themselves. Up to 20 first 
order equations can be input to the package at any one time. 
Note that at the end of each equation a semicolon must ap- 
pear except for the last equation which must be followed by 
two semicolons. The file EQIN is used as the input data 
file for the program. 

Variable names are restricted to 9 or less alphanumeric 
characters, the first of which must be alphabetic. Constant 
values are restricted to 20 or fewer characters. The fol- 
lowing variables names are restricted; III, JJJ, KKK EPS, 
IND, YO, YNEW, 2ZZB, ZZZZ1, ZZZZ2, 1 1 II, NIII, KKK1 , and 
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3-letter variable name beginning with T whose second letter 
is not A, and if a variable is used as a time derivative, 
for example X., then both X and DX are also reserved. The 
name T is reserved for the independent variable. 

The constant definition section consists of predefined 
user constants. The names constantl, ..., constantn are 
unique variable names, and valuel, ..., valuen are numbers. 
These numbers can be integers, simple real numbers, or 
exponential real numbers. Their form is identical to FOR- 
TRAN constants, for example: 1, 10.3, -.3, -1.976E-5, and 
.43795E-15 are all valid. Recall that the constant defini- 
tion becomes a data statement in the FORTRAN subroutines so 
that naming conventions for FORTRAN variables must be fol- 
lowed. 

The differential equation specification section con- 
sists of a series of first order ordinary differential equa- 
tions in the independent variable T, each equation of the 
form: 

variablel. = equationl;. (5) 

All variables used as derivatives must be real, i.e., begin 
with A. .H or O..Z. Only derivatives may appear on the 
left-hand side of the equations; derivatives may not appear 
explicitly on the right-hand side. The right-hand side may 
contain constants from the constant definition, numbers, 
variables defined by differentiation on the left hand side, 
operators and predefined user functions listed below, the 
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independent variable t, parenthesis, and the terminating 
semicolon. Note that the syntax of the right-hand side will 
be identical to FORTRAN assignment statements. 

The standard operators and 4 predefined user functions 
are available in this package. These operators are +, /, 
*, and **. Operands for these operators may be constants, 
expressions, the variable t, or series (variables defined by 
differentiation) with the following exceptions: 

1) t, when appearing by itself may be raised only to integer 
powers 

2) constants may not be explicitly divided by t or powers of 
t 

3) series may not be explicitly divided by t raised to 
powers greater than 1. 

The four predefined user functions are SIN, COS, LOG (natu- 
ral logarithm), and EXP. The arguments for these functions 
must be completely enclosed in parenthesis and may be con- 
stants, t, expressions or series variables with one excep- 
tion; LOG ( t ) or LOG(t**power) are not allowed. 

The reason for these restrictions is inherent in the 
method used for solution of the problem itself. In the 
first case LOG(t) does not have a power series defined about 
t=o so it along with LOG( t**power) are not allowed; however, 
LOG( 1+t ) or LOG(constant + t**power), e.g., are allowed. The 
reason that t may be explicitly raised only to integer 
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powers is that the program considers each term to be a con- 
stant or a power series. If, however the term t**l/2*SIN(t)**3/2 
appears one can simply rewrite it as( t*SIN(t) )**l/2*SIN(t) 
which is a legal expression. In most cases this problem can 
be avoided by rewriting the input. Since division by power 
series with zero constant term is not defined series may not 
be divided explicitly by t raised to a power. However, 
division by t by itself is defined so that x/t**2 could be 
rewritten as (x/t)/t, or x**3/2/t**l/2 as (x/t)**l/2**. 

Note that COS(t)/(l + t**2) or x/(8.7 + t) are no problem. 
Finally, the reason that constants may not be divided by t 
is that power series do not contain terms in inverse powers 
of t. To avoid this problem input can simply be rewritten, 
for example (const/t)*EXP(x) would become const* ( EXP (x)/t) 
or (3/t**2 )*COS(t) would become 3 . *(COS(t)/t )/t. 

Other than these restrictions on inputs stated above 
two other problems exist with respect to power series opera- 
tions. The first is the problem encountered with division 
by powers of t; a power series used as a divisor must have a 
nonzero constant term. The second problem occurs when 
raising a power series to a power other than 2 in which case 
the series must also have a nonzero constant term. These 
problems will show up when subroutine SOL is executed and a 
division by zero error message will be generated. 
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The above problems occur in only a very limited number 
of equations; these are not in the author's opinion major 
stumbling blocks for use of this part of the software pack- 
age. Although the user may be inconvenienced by having to 
rewrite some of his input due to these restrictions, the 
generation of SOL and DIFFUN by this program are still a 
great savings in terms of time and effort by the user. 



EXAMPLES OF INPUTS 


a) X' = -X-2Y+X**2*SIN(T) 

Y‘ = 5X -(T+3 )/(T+4)*Y 
X(0) = Y(0 ) = 1 

input would be 

X. = -X-2.*y+X**2*SIN(T); 

Y. = 5 . *X- ( T+3 )/ ( T+4 . ) *Y ; ; 

YO(l) 1., YO(2 ) = 1. 

b) X' = Y 

Y' = ( 5*EXP(T/TS ) - 1 - 3/2*Y**2 )/S + (A*Y+D)/X**2+C/X**(3*G+1) 
G = 1.4, A = 4.0E-2, TS = 29, D = .1456, 

C = 1+D 

X(0) = 1, Y(0) = 0 
input would be 

/* G = 1.4, A=4. 0E2 , TS=29 . , D=.1456, C=1.1456*/ 

X. =Y; 

Y . = ( 5 . *EXP ( -T/TS ) - 1. -(3./2. )*Y**2)/X - 

(A*Y+D)/X**2 + C/X**(3.*G+1. );; 

YO(l )=1 . , YO (2) = 0. 

c) D= . 349066T 
A=ATAN(W/U) 

THETA ' =Q 

U'=C*S IN (THETA) -W*Q + (Cl+C2*A)*(U**2+W**2 ) 

+ C3*Q*S&RT(U**2+W**2) 

W' = C*COS (THETA) + U*Q + (C4+C5*D+C6*A)*(U**2 
+W**2 ) + C7*Q*SQRT(U**2+W**2) 
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Q* = (U**2+W**2) * (C8+C9*(A+D) ) + (C10*W'/U 
+ C11*Q) *SQRT (U**2*W**2 ) 

C = 32.16, Cl = . 5018521E-9, Cl=0.1192125E-5, 

C3=.42110675E-3, C4=. 14598254E-3, C5=0.46345531E-5, 
C6=.45006065E-5 ; C7=-.69718949E-3, C8=.587215E-8, 
C9=.73872E-6, C10=0.3800645E-5, C11=.7422436E2, 

U(0)=660. 16187, W(0)=5. 74626, Q(0)=0, THETA(0)=5.3 . 

Notice that W' appears on the right-hand side, ATAN is not a 
user function, and A and D are not in the form of deriva- 
tives; however, we can substitute 
D.=. 349066, D(0)=0. 

A.=(U*W.W*U.)/(U**2+W**2), A(0)=ATAN(W(0)/U(0) ) 
the input becomes 

/* C=32 . 16, Cl=. 5018521E-9, C2=0.1192125E-5, 
C3=.42110675E-3, C4=.14598254E-3, C5=0.46345531E-5, 
C6=-.45006065E-5, C7=. 69718949E3 , C8=. 587215E8, 
C9=.73872E=6, C100.3800645E-5, Cll=.7422436E-2 */ 

U.=C*S IN (THETA) - W*Q + (C1+C2*A)*(U**2+W**2 ) 

+ C3*Q*(U**2+W**2 )**.5; 

W.=C*COS( THETA) + U*Q + (C4+C5*D+C6*A)*(U**2+ W**2) 

+ C7*Q*(U**2+W**2)**.5; 

Q.=(U**2+W**2)*(C4+C5*D+C6*A) + (C10*(C*COC(THETA) 

+ U*Q + (C4+C5*D+C6*A)*(U**2+W**2 ) 

+ C7*Q*(U**2+W**2) **. 5)/U+Cll*Q)*(U**2+W**2 )**.5; 
THETA. =Q; 
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A. = (U*(C*COS(THETA) + U*Q + C4+C5*D+C6*A) * (U**2 

W**2) + C7*Q* (U**2+W**2 ) ** . 5 ) - W*(C*SIN( THETA) 

- W*Q + (C1+C2*A)*(U**2+W**2) + C3*Q*(U**2+W**2 ) 

** • 5 ) )/(U**2+W**2 ) ? 

D. s. 349066; ; . 

Y0(l)=660. 18167, YO(2)=5.74626, YO(3)=0., 

YO(4)=5.3, Y0(5) S ATAN(W(0)/U(0) ), Y0(6)=0. 
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OUTPUTS 


There are two fortran subroutines output. The first is 
subroutine DIFFUN(T,Y,DY) which appears in source form on 
file DIFFUN. The variables Y and DY are 20 element arrays. 
The second is subroutine SOL(T, YO, YNEW, IND) which appears on 
file SOL. The variables YO and YNEW are also 20 element 
arrays. One important fact to notice is that the ordering 
of the input equations determines the coefficients in Y, DY, 
YO, YNEW, beginning at 1 corresponding to each equation, for 
example if: 

X.=Y; 

Y.--X; ; 

were entered X and DX would correspond to Y(l), DY(1), 
YO(l) , and YNEW(l); Y and DY would correspond to the sub- 
script 2 for both input and output values from the sub- 
routines . 

The purpose of subroutine SOL is to solve the system of 
equations analytically so that the first 10 significant 
digits are correct. The initial values at T=0 for the 
equations in order of occurrance are placed in YO and along 
with a value of T are input to SOL. The solution to the 
equations at this specific T are output in order of occur- 
rence in YNEW. The value of IND at output is zero if the 
radius of convergence of all series is (-1,1) and greater 
than zero if the radius of convergence of any series is less 
than (-1,1). Subroutine DIFFUN is explained in chapters 2 
and 4. 


61 



NOS CONTROL CARD EXAMPLES 


BATCH JOB FROM CARDS 
JOB. 

USER. 

CHARGE. 

ATTACH, SERIES/UN=USERNAM. 
COPYBR, INPUT, EQIN. 

REWIND, EQIN. 

SERIES , F-105000 . 

REPLACE, SOL. 

REPLACE , D I FFUN . 

*EOR 

DATA ON CARDS 



INTERACTIVE SUBMISSION AFTER CREATION OF FILE EQIN 


ATTACH, SERIES/UN=USERNAM 

GET, EQIN 

SERIES, F=105000 

REPLACE, SOL 

REPLACE, DIFFUN 

*EOF 

If running interactively using the SUBMIT command add the 
JOB, USER, and CHARGE cards to the above control card se- 
quence . 
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ERRORS 


When a compile time error occurs due to a mistake in 
the user input a dump is written by the program onto file 
OUTPUT. This dump contains the following information: 

a) a general statement such as ILLEGAL OPERATOR to give 
the user some idea of what type of error has occurred 

b) a listing of the input characters up to the point 
where the error was detected 

c) a dump of the symbol table that shows which variable 
and constants have already been b ered 

d) a listing of the codetable which shows how far into 
the equations the parser has progressed. 

The written messages in a) are neither completely descrip- 
tive nor always correct in the description given. If the 
parser becomes confused the detection of an error may occur 
well after the actual error has been made. If, for example, 
an opening parenthesis is missing in the expression the 
input may still remain legal until the extra closing paren- 
thesis is encountered. The best way to find a compile time 
error is to study the listing b) near the end of the char- 
acter listing. The symbol table and code table can also be 
helpful in error detection. When looking at the listing of 
the symbol table all the variable names and operators that 
have been encountered or generated will be listed. Next to 
each name is a number in parenthesis which tells the type of 
the variable: 
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0) user function 

1) vector or series variable 

2) constant 

3) operator or delimeter 

4) semicolon 

5) differentiated variable 

6) independent variable t. 

Variables of type 1 and 5 should be declared as arrays in 
subroutine SOL. The code table is listed as a series of 
quadruples describing the operations which are to be per- 
formed in order of occurrance. By cross-referencing between 
the code and symbol tables many logical errors will be 
easier to find. Note that. only the first 10 characters of 
numbers will be listed. Finally if the message RUNTIME 
STACK OVERFLOW appears in the user day file simply increase 
the F parameter on the SERIES control card. 
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LIMITED EXTENSIONS 


The biggest thing to remember is that this part of the 
package is only a compiler to generate the FORTRAN code. 
All the variables need not be defined in the input at com- 
pile time when the subroutines are generated. If one has 
any function which does not depend on the differentiated 
variables (VAR.) one could add it to subroutine SOL or 
DIFFUN after creation. This could be done using a statement 
function in the subroutines themselves or by use of a func- 
tion subprogram used in conjunction with SOL and/or DIFFUN. 
This extension is especially useful for functions of t. As 
a simple example let us redo the variable D in the previous 
example C. Recall that: 

D = . 349066T 

one could simple define the statement function, 

H(R) = . 349066*R 

and set D=H(T) at the beginning of the program. If any of 
the functions arguments are the dependent variables then one 
must define the functions as an initial value problem with 
initial conditions at t=0 such as A=ATAN(W/U) in example C. 
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4. USER'S INSTRUCTIONS FOR CLMP 
The user of the closed Loop Modelling package (CLMP) is 
faced with the following system of equations: 

Yi(t) = f i (t,y(t),u(t)) > i=l # 2,...,n (1) 

where u(t) is an input function of three possible types: 

(i) Damped sine wave 

(ii) Spike function 

(iii) Step function 

In order to use CLMP, the user must supply the fol- 
lowing information: 

(i) A system of equations 

(ii) An interval [t Q , t n ) on which y(t) is to be deter- 
mined 

(iii) Initial values for y(t Q ) 

(iv) Parameters of u(t) 

(v) A fixed-step-size integration routine for solving 
differential equations. 

A System of Equations 

The system of equations (1) will be stored on file 
DIFFUN by the user, as outlined in section 3. 

The system will be solved for y(t) for t in (t Q , t n ), 
where t Q and t must be specified by the user. It is also 
essential that the values of y(t) at the point t = be 
specified. 


67 




Should the user wish to write his own DIFFUN routine, 
rather than use the SERIES package for this, he may do so. 
The format to use is 

SUBROUTINE DIFFUN (T,Y,DY) 

DIMENSION Y(20) , DY(20) 

The routine simply computes DY(I) = f^ (t,y(t),u(t) ), 
as in (1). This routine should be compiled and stored on 
permanent file DIFFUN, as follows: 

FTN. 

SAVE , LGO=D I FFUN . 

7/8/9 

SUBROUTINE DIFFUN 

6/7/8/9 

Parameters of u(t) 

The first parameter is IUTYPE. If IUTYPE = 1, then 
u(t) is a damped sine wave with equation 

u(t) = e" at sin(pt+(j) ) 

If IUTYPE = 2, then u(t) is a spike function, and if 
IUTYPE = 3, u(t) is a step function (see figure 4.1). 

Other parameters which must be specified are a, b, 
UMIN, and UMAX, which simply indicate the domain and range 
of the function u(t). It should be noted that for IUTYPE = 2 
and IUTYPE = 3, the function u(t) is different from UMIN 
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only on the first one-fourth of the interval [a,b]. (i.e., 

[a,a*(b-a)/4] ) . 

Integration Routine 

CLMP provides the user with 7 standard fixed-step-size 
integration routines to choose from. They are: 

(1) Euler's method 

(2) Improved Euler (2-pt. Runge-Kutta) 

(3) 2-pt. Adams -Bashforth 

(4) 3-pt. Runge-Kutta 

(5) 4-pt. Runge-Kutta 

( 6 ) Runge-Kutta-Merson 

(7) 4-pt. Adams -Moulton. 

It should be noted that all of these methods are self- 
starting, except for (3) and (7), which use Improved Euler 
and 4-pt. Runge-Kutta to start, respectively. 

Should the user desire to test his own integration 
routine, instead of one of the above, he may do so. Such a 
routine, however, must be called NEXTPT and have the calling 
sequence : 

SUBROUTINE NEXTPT(T,Y,N, STEP,KTR) 

NEXTPT, when given the values for y(t), should compute 
y(t+h), and return this value in y. Here, h is a fixed- 
step-size specified by the variable STEP. CLMP requires 
that STEP < (t n - t Q )/32. The array Y is dimensioned as 
Y(20 ) , and y^(t) is contained m Y(I). The variable N 
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" 'T 




denotes the number of equations in the system, and cannot 
exceed 20. The variable KTR may be used for whatever pur- 
pose is necessary. When the first call to NEXTPT is made 
from CLMP, the value KTR=0 is transmitted and remains at 
KTR=0 until changed in NEXTPT. No other information is 
passed through KTR. 

A sample problem and terminal session follow: 

Sample Problem 

The following equations govern the longitudinal motion 
of a jet fighter [Steinmetz et al. ] : 


a = -g sin e-wg ♦ ^ [<C x ) #t(4j + 


, 2 „- 


a = g cos 9 + uq + + (C^fS-Sj)] 


4 = [(c n > « T ,« T 


+ C (a-a m ) + C ££ + C 
m ' T' m. 2V nr 
a a q 


+ C, 


m fi * 6 “ 6 T^ 


0 = q 

a_ = w - g cos 0 - qu 

JU 

a ~ H 

a = tan-1 (~) 
v = 
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where all C values (C , etc.) are treated as constants. 

“a 

These equations were set up in DIFFUN as shown in 
Figure 4.2. The input for 6 was taken from the function 
u(t) with parameters IUTYPE, UMIN, UMAX, A, and B read as 
data. The data for the test run was as follows: 


N 

TO,TN 

(YO(I), 1=1, N) 
IUTYPE 

UMIN, UMAX, A, B 


4 

0 ., 10 . 

660.18167, 5.74626, .5, .09251 
2 

0., 0.174533, 0., 4. 


STEP 


.1 


Here, Y(l) is y, Y(2) is w, Y(3) is Q, and Y(4) is 0. 

The values for C w , C_ , etc . , were extracted from data 

x a a 

given in [Steinmetz, et al.], and are shown in DATA state- 
ments in Figure 4.2. 
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Figure 4.2 Input differential system for CU(P. 


SUBROUTINE DIFFUN(T, Y,DY) REAL Y(13,20), CY(20) 

C *********************************************************** 
C* THIS REPRESENTS THE LONGITUDINAL STABILITY EQUATIONS 
C* FOR A JET FIGHTER. THE INPUT FOR STABILATOR DEFLECTION 
C* IS FROM THE FUNCTION U(T), WITH PARAMETERS: IUTYPE=1, 

C* UMIN=- . 05 , UMAX=.065, A=0. , B=10. 

£*********************************************************** 

REAL M, IY. 

DATA S, G, M, RHO, CBAR, AT, DT, IY/530., 32.158, 

1285.5, 2.125E-3, +16.04, 0., 0., 129609/ 

DATA CXATDT, OXA, OXQ / . 0158, 0890, -3.92 / 

DATA CZATDT, C2A, C2Q, CZD / -.121, -3.36, -6.49, .346 / 

DATA CMATDT, OMA, CMADOT, CMQ, CMD / -.00117, -.106, 

-1.66, -1.66, -1.66, .576/ 


UU=Y (1,1) 

W=Y (1,2) 

Q=Y(1,3) 

THETA=Y (1,4) 

V=SQRT(UU*UU+W*W) 

ALFA=ATAN ( W/UU ) 

DELTA=U ( T ) 

Cl=( .5*Q*0BAR)/V 
C2= . 5*RHO*V*V*S 
A=ALFA-AT 
D=DELTA-DT 

DY ( 1 ) =-G*S IN ( THETA ) -W*Q+ ( 02/M ) * ( 0XATDAT+0XA*A+0XQ*01 )DY( 2 ) 
=G*C0S ( THETA ) +UU*Q+ ( C2/M ) * ( CZATDT+CZA*A+CZQ*C1+CZD*D ) 

ALFD0T=DY ( 2 )/UU 

AZ=DY ( 2 ) -G*C0S ( THETA ) -Q*UU 

DY ( 3 ) = ( C2*CEAR/1Y ) * ( OMATDT=CMA*A 
+ =( . 5*CMADOT*ALFDOT*CBAR)/V+CMQ*Cl+CMD*D) 

DY(4) = Q 

RETURN END 
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Running a Job on CLMP 

As mentioned earlier, the user may either supply the 
system of equations (1) directly (sect. 3) or as SUBROUTINE 
DIFFUN. See Figure 4.3. Procedures for running CIMP on a 
CYBER 172 computer with NOS 1 operating system are given as 
follows: 

Equations entered directly: 

(Jobcards) 

COPYBF( INPUT, TAPE6) 

(CLMP commands) 

6/7/8/9 

Equations entered in DIFFUN: 

( Jobcarbs) 

FTN,B=DIFOBJ. 

7/8/9 
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FIGURE 4.3 SAMPLE TERMINAL INPUT 


WHAT IS THE RATE OF YOUR TERMINAL IN CHARACTERS PER SECOND? 
>>30 

HOW MANY EQUATIONS ARE IN THE SYSTEM? 

?4 

GIVE INITIAL AND FINAL VALUES OF T 
? 0 ., 10 . 

GIVE INITIAL VALUES FOR 
Y(l) 

? 660.18167 
Y<2) 

? 5.74626 
Y(3) 

? .5 

Y(4) 

? 5.3 WHAT TYPE OF NOISE WILL BE INPUT (U(T) )? 

ENTER 1. FOR DAMPED SINE WAVE 

2. FOR SPIKE FUNCTION 

3. FOR STEP FUNCTION 

? 1 

GIVE MINIMUM AND MAXIMUM VALUES FOR U(T) 

? -.05,. 065 

WHAT INTERVAL WILL U(T) BE DEFINED FOR? 

? 0., 10. WHAT STEP SIZE WILL BE USED FOR T? 

?.l 

WHICH INTEGRATION ROUTINE DO YOU WISH TO COMPARE WITH DIFSUB? 
ENTER 1 FOR EULER 

2 IMPROVED EULER 

3 2-PT ADAMS -BASHFORTH 

4 3-PT RUNGE-KUTTA 

5 4-PT RUNGE-KUTTA 

6 RUNGE-KUTTA-MERSON 

7 4-PT ADAMS -MOULTON 

?6 
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SUBROUTINE DIFFUN 


END 


7/8/9 

(CLMP commands ) 

6/7/8/9 

The section described simply as (CLMP commands) now 
follows : 

GET, A=CLMP 
UVAUB, PLOTIO. 

ENTER. /LOAD( A, DIFOBJ )/NOGO, ABS . 

LABS , ABS . 

( options ) 

Since values are output on TAPE6 and TAPE7, the user 
may desire to: 

(i) Save TAPE6 and TAPE7 for further use 

(ii) List TAPE6 and/or TAPE7 at the line printer. 

(iii) Have the values on TAPE7 punched at the card punch 
These may be accomplished by entering the following commands 
in the section referred to as (options): 

(i) SAVE , TAPE6 , TAPE7 . 

(ii) f .EWIND , TAPE6 , TAPE7 . 

COPYBF , TAPE6 , OUTPUT . 

COPYSBY , TAPE7 , OUTPUT . 
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(iii) REWIND, TAPE7. 

C0PYBF,TAPE7, PUNCH. 

Notes 

The Fourier coefficients given are of the form: 

f<t) = ^ + k=l Cr C0S <d,t + ^o> } f ° r t£lt °' (2) 

Only the first N coefficients are computed however, 
where N is a suitable value between 32 and 64, chosen by 
CLMP. They are stored in an array in the following fashion: 

ARRAY* 1) = c q /2 

ARRAY* 2) = 0. 

ARRAY* 3 ) = c x 

ARRAY * 4 ) = d x 
• 

* 

ARRAY (2N-1 ) = C n _ ] _ 

ARRAY (2N) - d n _ 1 

Using the construction in *2) gives not f(t), however, 
but instead, f(t+t Q ). Hence, in order to compute f(t) for t 
in [t Q , t n J, one must first compute t* = t - t Q , and then 
compute 

N 2knt* 

f < t) = ARRAY* 1) + I ARRAY * 2 k-1 )*COS* ARRAY* 2k) + , r- -r->) 

i=2 

Smoothing 

The FORTRAN subroutine RFFT is used by CLMP to compute 
Fourier coefficients for f*t). RFFT requires 64 points to 


77 



be entered as tables values, whereas CIMP will supply only N 
values, with 32 < N < 64. The remaining 64- N points are 
assigned value zero. The Fourier coefficients computed by 
RFFT yield a function f*(t) with the following property: 
Given f(t) and the interval [t Q , t n ),f(t^) - f*(t^) * 
constant for 1 < i < and f(t^) - f*(t^) has alternating 
signs for consecutive values of i. The value of the con- 
stant is unknown, but, because of the alternating signs, can 
be virtually eliminated graphically as follows: 

Given t^ and f (t^), 1 < i < n, 

compute T i = + ^i+i)/ 2 anci 

(Fftj) = { f * ( ) + £*(t i+1 ))/2. 

Notice that f(t^) * (f(t^) + f(t^ +1 ))/2 as well. Then by 
plotting F(t^) at t^, for 1 < i < N - 1, a smooth graph can 
be drawn that closely approximates f(t). 

By way of an example, examine figures 4.4, 4.5 and 4.6. 
Figure 4.5 shows the graphs of Y(l) computed from two 
separate sets of Fourier coefficients, before sm*. hing. 
Notice the undesirable effects of oscillation. Figure 4.6 
shows the same graphs, after smoothing. Notice how closely 
the DIFSUB and NEXTPT graphs coincide on the smooth portion 
of the graph (i.e., away from discontinuities in the first 
derivative). This gives a very good view of what Y(l) looks 
like on that interval. 
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APPENDIX A 
NUMERICAL METHODS 


ABK Predictors 

K=1 y„ • Vi + h y'n-i 

K=2 y n = Vi ♦ h / 2 ( 3 y'n-i - y 'n-2> 

K=3 y n = y n -l + h/12(23y- n . 1 -16y' n _ 2 + 5y' n . 3 ) 

K=4 y n = y Dml + 11/24(557'^ -59y' n _2 + 37y' n . 3 -9y’ n _4> 

K=5 y n = - v n-l + V720(1901y' n . 1 -2774y' n _ 2 + 2616y' n _ 3 

-1274y' n . 4 + 251y' n _ 5 ) 

K=6 y n = y n . 1 + h/1440 (4277y' n . 1 -7923y' n _ 2 +9982y' n _ 3 
-7298y' n . 4 + 2877y' n _ 5 -457y' n _ 6 ) 

AMK - Corrector (used with ABK Predictor) 

K=1 y n = y n .j + h/2(y’ n + y' n _l> 

K= 2 y n = y n .i + h/12(5y' n + Sy^.2 -y’ n _2> 

K=3 y„ = y n-l + V24(9y' n + 19y’ n . 3 -Sy' n _2 + y'n-3* 
K=4 y n = y n -l + h/720(251y' n + 646y' n _ 1 -264y' n _ 2 + 
I06y' n . 3 -I9y' n _ 4 ) 

K-5 y n = y n . 3 + h/1440(475y' n + 1427y’ n _ 1 -798y' n . 2 + 
482y' n _ 3 -173y' n . 4 + 27y' n _ 5 ) 

BDF - Corrector (used with ABK Predictor) 

K=1 y n = y n -i + hy 'n 

K=2 3y n = 4y n _ 3 -y n . 2 «hy' n 

K=3 lly n = 18y n . 1 -9y n . 2 + 2y n _ 3 +6y' n 

K = 4 25y n = 48y n _ 1 -36y n _ 2 + 16y n . 3 -3y n . 4 +12hy' n 
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75y n-4 

■ 225y„. 4 ♦ 72y n . S 


RK2 - one step 

*o * h£ <y n -i> 

K i = hf( y n -i + v 2 ) 

^ =y n -i + K i 

RK4 - one step 

K 0 * hf <y„-i> 

h = hf (y n -i + v 2 > 
k 2 = hf <y n -i + v 2 > 

K 3 = hf <y n -i + k 2 > 

y n = y n-l + h /6 < K 0 + 2K 1 + 2K 2 + K 3 i 


K=5 137y n = 300y n . 1 - 300y n _ 2 + 200y n _ 3 - 
+ 12y n-5 + 60hy ’i; 

K=6 147y n = aeoy^ - 450y n _ 2 + 400y n _ 3 
' 10y n-6 + “hy'n 
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APPENDIX B 


THEOREM PROOFS 

Page 

Theorem 1 B1 

Theorem 2 B2 

Existence Theorems B3 
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Definition. A linear k-step formula satisfies 

W 0 * £ fy*-/ + /• /„_,)] . 

Definition. A one-leg k-step formula corresponding to (9) satisfies 

(10) 0 - tv.-, + i 7 t v.-/) • 

where s = o(l). Without loss of generality, set s = 1. 

Theorem 1 (4). Let Y n be a sequence which satisfies (10), and let Y n = 

{?„} be such that 

(id y n - £ V-/ ■ 

/= o 

w/tere £ denotes the back shifting operator. Then Y n satisfies (9). Conversely, if 
Y n satisfies (9), then there exists a sequence Y n such that y n = o{E)y n , and y„ 
satisfies (10). 

Proof. Without loss of generality assume the system of equations is autonomous. 

Write (10) as p(£)>’ fl = - hf(o(E)y n ), n = 0, 1, 2 This together with (1 1) 

implies p(E)y n = p{E)o(E)y n = - ho(E)f(y „), which implies that y n satisfies (9). 

For the converse, Euclid’s theorem on polynomials with no common divisors 
implies the existence of polynomials P, Q which satisfy P(x)o(x) + Q(x)p(x) = x m , 

0 <m <k. Writing (9) as p(E)y n = - ho(E)fiy n ) and setting v„ = 
E~ m (P{E)y n - hmf(S' n )) gives 

o{E)y n = E~ m {P{E)o{E)y n + Q (E)p{E)f(y n )) = E~ m (P(E)o(E) + Q(E)p{E))y n , 
which gives o(E)y n = y n . Next, set 

p(E)y„ = - E~ m h(P(E)o(E) + Q(E)p(E))f(y n ) = - hf(y n ) = - hf(o(E)y n ), n > m, 

and Y n satisfies (10), proving the converse. 

This shows that Y n given by the one-leg fc-step formula will have similar 
stability properties to its corresponding linear A-step sequence y. Dahlquist [4] lias 
described a discrete Liapunov function V G ,t,h which, applied to a sequence Y n , 
characterizes the stability of that sequence generated by a nonlinear system (1). Let 
v a,i,hi Y n) = £/=/S/=i«// >’«_/+!>. where G is a positive definite, l x l 

matri... The structure of G assures that V G t h is positive for Y n # {0}. 


B1 



Theorem 2. If V cth (Y n ) * c for the symmetric positive definite matrix G, 
then there exists a symmetric, positive definite matrix G, dependent only on G and 
o(x), such that Vg IH (Y„) - c, where y n ■ <;(£>„ are the elements of Y n and Y n . 

Proof. Without loss of generality, consider a system of only one equation / * 
f(y, t) generating the sequence Y n = {y„_* +) , • * • >?„)■ Since y H andy„ are 
related only by the k + 1 coefficients of o(x), replace the sequences by the vectors 

= O',. y»- k + 1 )* “d w n = (y n > 1 >• where * ** 

the transpose operator. Let & be a (k + 0 by (Ar + /) matrix consisting of G in the 
upper l by l partition, and 0 elsewhere. Then I'Ov,,) - w*Gw„ = c > 0, and 
Y(w’„) - w'*G'w' n - c. 

Define S such that w H = Sw' n , thus 


ho 

f 0 b 0 


S * 


b k 0 ••• 0 ^ 
**-■ *, • ■ • 0 



h o *i *./ 


is an / by / + k matrix. Then since G’ is of rank l because G is positive definite, 
there exists a singular value decomposition of (f = UFV* for U, V l + k by / + k 
unitary matrices, and F = (f q ) for D an / by / diagonal matrix of singular values of 
G. Thus, there exists an / by / matrix G such that G' = S*GS. This is seen by letting 
S have the singular value decomposition £/,(S |0) V f , where 2 is an / by / diagonal 
matrix. Then 


d = u,( 2 -* \q)v*ufv*v *&- 1 \o)*u*. 

If b 0 , b v . . . , b t _ , are all zero, a similar argument can be made using an / 
by / + k - i matrix S where b t is the coefficient of lowest index / such that b t & 0, 
and 



Thus, there exists a G dependent only on o(x) such that Yg , h (Y n ) = c for all c, 
which was to be shown. 
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Theorem 4. If the equilibrium p(t) of f(y, t) has a Liapunov function of the form p(y, t) - y*Qy 
for a positive definite matrix Q, y* the transpose of y, and f(y, t) is continuously differentiable 
on a convex domain D whose boundary 3D is (Mined by v'(y 0 , t«) ■ 0, and if f(y, t) has the 
property on A that (y, - yi)*Q(f(y„ t) - /(y 2 , /)) s - yj|<? 2 where m *0 and x*Qx -Mo*, 
then for any point ft - 9(t t ) in the interim' of A, an (/, /, h)-stability region can be constructed 
using rays from ft, provided h s h<fl,f). 

Proof. Note that the solution from r 0 to r«+h is spiraling in from 3D toward the 
equilibrium, which is inside a circle lyfO-yJjSftAOW where A/«max/'(r), since the 

ttD 

hypotheses and (8) gives |yft) - y(//)| * exp (jMfbUti - Hhli By definition, AV, M (Y,)~ 
which occurs when 0 - v(y h /,) - v(yo, f 0 ) - (y t - yo)v'(z) by the mean 
value theorem. Since c'(z) m Q on the boundary 3D, along any ray from ft one can find y(ti) 
generated by a y 0 on 3D, provided hi < hoi is small enough that the trajectory from ft> to ft is 
entirely in A 

The boundary of (/, /, A)-domain of attraction A mm be constructed of all such rays. The 

(/, /, /^stability region A has a boundary of all points such that V, JA « e* ■ min V Uk , which 

•o 

can also be constructed. 

Theorem 5. If the hypotheses of theorem 4 are met and (3l3y)(i l> f(z)ldl p ) is continuous in A 
then for any ft in the interior of A, an (/, /, h) stability region can be constructed for a pth order 
one-leg Jt-step method, for Hsh^l,f). 

Proof. A V « |y,| 0 2 - |y(f, - lh%<} where 

yi ~ offd,.,) + hf(^ bjy(t H ), ^ bjt ,., ) 

and y(ti) m yi + K p h'* t f <,) (z), the truncation error formula. If the truncation error varies 
continuously as y(t o) varies along 3D, Then two solutions y» and z/ generated by yo, r« on 3D 
have the property that y ( -*$ uniformly as yo-*r», and a smooth curve 3D" exists on the 
boundary of the (/, /, h)-domain of attraction. Similar arguments show the existence of the 
(/, /, hl-stability region. 
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